library(tidySummarizedExperiment)
library(GEOquery)
library(tidyverse)
library(ggpubr)

jzl_list <- read_csv('Archive/FPP_TRPV3/data/JZL_list.txt')

fpp_gene <- read_csv('Archive/BD-MassSpec/data/fpp_genes.txt') |>
  pull(1)

care_list <- c(jzl_list[[1]], jzl_list[[2]], fpp_gene)

# GEO2R top table ----------
peng_zhu_2021 <- read_tsv('Archive/FPP_TRPV3/data/GSE172423.top.table.tsv')

peng_corey_2009 <- read_tsv('Archive/FPP_TRPV3/data/GSE18527.top.table.tsv') |>
  distinct(Gene.symbol, .keep_all = TRUE)

sdeg2021 <- peng_zhu_2021 |>
  filter(ID %in% care_list & adj.P.Val < 0.05)

peng_zhu_2021 |>
  mutate(DEG = case_when(
    logFC > 0 & adj.P.Val < .05 ~ 'Downregulated in lesion',
    logFC < 0 & adj.P.Val < .05 ~ 'Upregulated in lesion',
    .default = 'Non-significant'
  )) |>
  ggplot(aes(-logFC, -log(P.Value), color = DEG)) +
  geom_point() +
  geom_point(data = sdeg2021, color = 'red') +
  scale_color_manual(values = c('blue','grey','orange')) +
  theme_pubr() +
  labs_pubr() +
  theme(text = element_text(size = 18)) +
  ggrepel::geom_label_repel(data = sdeg2021,
                            aes(-logFC, -log(P.Value), label = ID),
                            color = 'red',
                            size = 5,
                            fontface = 'bold')

sdeg2009 <- peng_corey_2009 |>
  filter(Gene.symbol %in% care_list & adj.P.Val < 0.05)

peng_corey_2009 |>
  mutate(DEG = case_when(
    logFC > 0 & adj.P.Val < .05 ~ 'Downregulated in lesion',
    logFC < 0 & adj.P.Val < .05 ~ 'Upregulated in lesion',
    .default = 'Non-significant'
  )) |>
  ggplot(aes(-logFC, -log(P.Value), color = DEG)) +
  geom_point() +
  geom_point(data = sdeg2009, color = 'red') +
  scale_color_manual(values = c('blue','grey','orange')) +
  theme_pubr() +
  labs_pubr() +
  theme(text = element_text(size = 18)) +
  ggrepel::geom_label_repel(data = sdeg2009,
                           aes(-logFC, -log(P.Value), label = Gene.symbol),
                           color = 'red',
                           size = 5,
                           fontface = 'bold')

# load in GEO set -----------
tseobj <- getGEO("GSE172423",
               filename = 'Archive/FPP_TRPV3/data/GSE172423_series_matrix.txt.gz',
               GSEMatrix = TRUE, 
               AnnotGPL = FALSE) |>
  makeSummarizedExperimentFromExpressionSet()

fpp_expr <- tseobj |>
  select(c(.feature, title, exprs, time.point.ch1)) |>
  filter(exprs > 0 & .feature %in% sdeg2021$ID) |>
  mutate(log_exprs = log2(exprs)) 

p_coord <- fpp_expr |>
  group_by(.feature) |>
  summarise(y = max(log_exprs) +0.2) |>
  pull(y)

sdeg2021_prism <- sdeg2021 |>
  arrange(ID) |>
  mutate(adj.P.Val = signif(adj.P.Val, 2), .feature = ID) |>
  add_column(group1 = 'ctrl',
             group2 = 'lesion',
             y.position = p_coord)

fpp_expr |>
  mutate(time.point.ch1 = fct_relevel(time.point.ch1, 'ctrl', after = Inf)) |>
  ggplot(aes(time.point.ch1, log_exprs, fill = time.point.ch1)) +
  geom_boxplot() +
  facet_wrap(~.feature, scales = 'free') +
  labs(y = 'log-expression',
       x = 'time post-infection',
       fill = 'time post-infection',
       title = 'Human genital skin biopsy after HSV-2 infection') +
  theme_pubr() +
  labs_pubr()+
  stat_pvalue_manual(sdeg2021_prism, label = 'adj.P.Val', inherit.aes = FALSE) +
  scale_fill_manual(values = c('orange','yellow','red','blue'))+
  scale_y_continuous(expand = expansion(mult = c(0.1,0.2)))

# 2009 data boxplot -----------
tseobj <- getGEO("GSE18527",
                 filename = 'Archive/FPP_TRPV3/data/GSE18527_series_matrix.txt.gz',
                 GSEMatrix = TRUE, 
                 AnnotGPL = FALSE) |>
  makeSummarizedExperimentFromExpressionSet()

probe_anno <- peng_corey_2009 |>
  select(c(ID, Gene.symbol))

fpp_expr <- tseobj |>
  select(c(.feature, title, exprs)) |>
  right_join(probe_anno, by = c('.feature' = 'ID')) |>
  filter(exprs > 0 & Gene.symbol %in% care_list) |>
  filter(str_detect(title, 'control|lesion') & !str_detect(title, 'for')) |>
  mutate(log_exprs = log2(exprs),
         group = str_extract(title, 'control|lesion')) |>
  filter(Gene.symbol %in% sdeg2009$Gene.symbol)

p_coord <- fpp_expr |>
  group_by(Gene.symbol) |>
  summarise(y = max(log_exprs) +0.2) |>
  pull(y)

sdeg2009_prism <- sdeg2009 |>
  arrange(Gene.symbol) |>
  mutate(adj.P.Val = signif(adj.P.Val, 2)) |>
  add_column(group1 = 'control',
             group2 = 'lesion',
             y.position = p_coord)

fpp_expr |>
  ggplot(aes(group, log_exprs, fill = group)) +
  geom_boxplot() +
  labs(y = 'log-expression',
       title = 'Human skin biopsies during the mucocutaneous HSV-2 reactivation') +
  theme_pubr() +
  labs_pubr()+
  stat_pvalue_manual(sdeg2009_prism, label = 'adj.P.Val', inherit.aes = FALSE) +
  facet_wrap(~Gene.symbol, scales = 'free') +
  scale_fill_manual(values = c('blue', 'red'))+
  scale_y_continuous(expand = expansion(mult = c(0.1,0.2)))
